Section: Research Program

Modeling

The main activities of Magique-3D in modeling are the derivation and the analysis of models that are based on mathematical physics and are suggested by geophysical problems. In particular, Magique-3D considers equations of interest for the oil industry and focuses on the development and the analysis of numerical models which are well-adapted to solve quickly and accurately problems set in very large or unbounded domains as it is generally the case in geophysics.

• Explicit High-Order Time Schemes. Using the full wave equation for migration implies very high computational burdens, in order to get high resolution images. Indeed, to improve the accuracy of the numerical solution, one must considerably reduce the space step, which is the distance between two points of the mesh representing the computational domain. Another solution consists in using high-order finite element methods, which are very accurate even with coarse meshes. However, to take fully advantage of the high-order space discretization, one has to develop also high-order time schemes. The most popular ones for geophysical applications are the modified equation scheme  [85] ,  [100] and the ADER scheme  [91] . Both rely on the same principle, which consists in applying a Taylor expansion in time to the solution of the wave equation. Then, the high-order derivatives with respect to the time are replaced by high order space operators, using the wave equation. Finally, auxiliary variables are introduced in order to transform the differential equation involving high-order operators into a system of differential equation with low order operators. The advantage of this technique is that it leads to explicit time schemes, which avoids the solution of huge linear systems. The counterpart is that the schemes are only conditionally stable, which means that the time step is constrained by a CFL (Courant-Friedrichs-Levy) condition. The CFL number defines an upper bound for the time step in such a way that the smaller the space step is, the higher the numbers of iterations will be. Magique-3D is working on the construction and the analysis of new explicit time schemes which have either larger CFL numbers or local CFL numbers. By this way, the computational costs can be reduced without hampering the accuracy of the numerical solution.

• Implicit High-Order Time Schemes. Solving wave propagation problems in realistic media and in time domain is still a challenge. Implicit numerical schemes are nowadays considered as too expensive because they require the inversion of a linear system at each time step, contrary to the explicit schemes. However, explicit schemes are stable only when conditions on the discretization parameters are fulfilled, which can be very difficult to satisfy in realistic contexts and lead to very expensive simulations. These conditions become less dramatic or even disappear in some cases when using implicit schemes. Our goal is to construct, justify and optimize analytically original implicit schemes that seem accurate to solve specific difficulties coming from realistic problems. Several directions could be followed. First, we will continue to develop a methodology to construct high order implicit schemes for simple domains (conservative and homogeneous). For now (in [23] ) we have used the modified equation technique on the classical $\theta$-scheme, which leads to a parametrized family of numerical schemes that do not possess the same consistency error. Then, instead of choosing a time step that leads to a good precision for a given numerical scheme and spatial discretization, we reverse this standard reasoning and choose the best stable scheme, in the family of schemes that we just built, for a spatial and temporal given discretization. Stability is shown by energy techniques. It would be possible to continue this approach, leading to higher order schemes and better mastering the methodology. Crucial improvements to this work will be to adapt the methodology to dissipative media, heterogeneous media, realistic boundary conditions and model coupling. For instance, we aim at developing locally implicit schemes, for which the degree of implicitness would depend on the local characteristics of the media. Implicit-explicit schemes would be an application case of these new schemes, that could be used to optimize the global cost of simulation. Since computational efficiency is a priority, this theoretical seek will systematically be completed by the study of associated algorithms and their implementation on parallel architectures. We believe that locally implicit schemes will be well suited to the use of parallel algorithms.

• Asymptotic methods for ultra short laser pulses propagation In the long term goal of modeling an entire ultrashort laser chain, our first objective is to model the propagation of an ultrashort laser pulse in an isotropic third order nonlinear dispersive medium (as silica which is the material used for optical fibers or lenses). In other words, the optical index of the medium depends on the wavelength ( dispersion phenomenon) but also on the electromagnetic field's intensity in a cubic way (Kerr effect). A first intuition is to use Maxwell's equations coupled with additional equations for the optical index. Current computing facilities allow us to solve such equations in parallel on small domains and during short time intervals, for instance using MONTJOIE software. The use of asymptotical methods that take advantage of the pulse's brevity leads to a family of equations written as evolution equations in the propagation direction (among which the nonlinear Schrödinger equation), and solved in frequency domain, which are much easier to solve. However, ultrashort pulses have large spectra, which contradicts another hypothesis currently done in usual asymptotic methods. This is why new models have to be derived, as well as numerical methods to solve them.

  In fiber optics, the laser pulse propagates inside a waveguide called “optical fiber”, in which the transversal spatial repartition of the electromagnetic field can be shown to be a linear combination of eigenmodes. A first idea will be to generalize the results obtained in 1d (see 6.2.12 ) to this more realistic application. We have good reasons to believe that a very efficient model will be derived and will compare very well with the global Maxwell system. An ultimate validation will be obtained by comparing the numerical results with experimental data.

  Following this step, and in collaboration with CEA-CESTA, we wish to derive this kind of asymptotic models and associated numerical methods for general 3D open laser propagation.

• Finite Element Methods for the time-harmonic wave equation. As an alternative to Time-Domain Seismic Imaging, geophysicists are more and more interested by Time-Harmonic Seismic Imaging. The drawback of Time Domain Seismic Imaging is that it requires either to store the solution at each time step of the computation, or to perform many solutions to the wave equation. The advantage of Time Harmonic problems is that the solutions can be computed independently for each frequency and the images are produced with only two computations of the wave equation and without storing the solution. The counterpart is that one has to solve a huge linear system, which can not be achieved today when considering realistic 3D elastic media, even with the tremendous progress of Scientific Computing. Discontinuous Galerkin Methods (DGM), which are well-suited for $hp$-adaptivity, allow for the use of coarser meshes without hampering the accuracy of the solution. We are confident that these methods will help us to reduce the size of the linear system to be solved, but they still have to be improved in order to tackle realistic 3D problems. However, there exists many different DGMs, and the choice of the most appropriate one for geophysical applications is still not obvious. Our objectives are a) to propose a benchmark in order to test the performances of DGMs for seismic applications and b) to improve the most performant DGMs in order to be able to tackle realistic applications. To these aims, we propose to work in the following directions :

  1. To implement a 2D and 3D solver for time harmonic acoustic and elastodynamic wave equation, based on the Interior Penalty Discontinuous Galerkin Method(IPDGM). The implementation of this solver has started few years ago (see Section  5.1 ) for solving Inverse Scattering Problems and the results we obtained in 2D let us presage that IPDGM will be well-adapted for geophysical problems.

  2. To develop a new hybridizable DG (HDG)  [84] for 2D and 3D elastodynamic equation. Instead of solving a linear system involving the degrees of freedom of all volumetric cells of the mesh, the principle of HDG consists in introducing a Lagrange multiplier representing the trace of the numerical solution on each face of the mesh. Hence, it reduces the number of unknowns of the global linear system and the volumetric solution is recovered thanks to a local computation on each element.

  3. To develop upscaling methods for very heterogeneous media. When the heterogenities are too small compared with the wavelengths of the waves, it is necessary to use such techniques, which are able to reproduce fine scale effects with computations on coarse meshes only.

  We also intend to consider finite elements methods where the basis functions are not polynomials, but solutions to the time-harmonic wave equations. We have already developed a numerical method based on plane wave basis functions  [89] . The numerical results we have obtained on academic test cases showed that the proposed method is not only more stable than the DGM, but also exhibits a better level of accuracy. These results were obtained by choosing the same plane waves for the basis functions of every element of the mesh. We are now considering a new methodology allowing for the optimization of the angle of incidence of the plane waves at the element level.

  Last, we are developing an original numerical methods where the basis functions are fundamental solutions to the Helmholtz equation, such as Bessel or Hankel functions. Moreover, each basis function is not defined element by element but on the whole domain. This allows for reducing the volumetric variational formulation to a surfacic variational formulation.

• Boundary conditions. The construction of efficient absorbing boundary conditions (ABC) is very important for solving wave equations. Indeed, wave problems are generally set in unbounded or very large domains and simulation requires to limit the computational domain by introducing an external boundary, the so-called absorbing boundary. This topic has been a very active research topic during the past twenty years and despite that, efficient ABCs are have still to be designed. Classical conditions are constructed to absorb propagating waves and Magique-3D is investigating the way of improving existing ABCs by introducing the modelling of evanescent and glancing waves. For that purpose, we consider the micro-local derivation of the Dirichlet-to-Neumann operator. The interest of our approach is that the derivation does not depend on the geometry of the absorbing surface.

  ABCs have been given up when Perfectly Matched Layers (PML) have been designed. PMLs have opened a large number of research directions and they are probably the most routinely used methods for modelling unbounded domains in geophysics. But in some cases, they turn out to be unstable. This is the case for some elastic media. We are thus considering the development of absorbing boundary conditions for elastodynamic media and in particular for Tilted Transverse Isotropic media, which are of high interest for geophysical applications.

• Asymptotic modeling.

  During the last 30 years, mathematicians have developed and justify approximate models with multiscale asymptotic analysis to deal with problems involving singularly perturbed geometry or problems with coefficients of different magnitude.

  Numerically, all these approximate models are of interest since they allow to mesh the computational domain without taking into account the small characteristic lengths. this techniques lead to a reduction of the computation burden. Unfortunately, these methods do not have penetrated the numerical community since most of the results have been obtained for the two dimensional Laplacian.

  The research activity of Magique 3D aims in extending this theory to three-dimensional challenging problems involving wave propagation phenomena. We address time harmonic and time dependent problems for acoustic waves, electromagnetic waves and elastodynamic wave which is a very important topic for industry. Moreover, it remains numerous open questions in the underlying mathematical problems.

  Another important issue is the modeling of boundary layers which are not governed by the same model than the rest of the computational domain. It is rather challenging to derive and to justify some matching condition between the boundary layer and the rest of the physical domain for such multiphysical problems.

  More precisely, we have worked in 2013 on the following topics:

  • Eddy current modeling in the context of electrothermic applications for the design of electromagnetic devices, in collaboration with laboratories Ampère, Laplace, Inria Team MC2, IRMAR, and F.R.S.-FNRS;

  • Multiphysic asymptotic modeling of multi perforate plates in turbo reactors in collaboration with Onera.

  • Modeling of small heterogeneities for the three dimensional time domain wave equation. This reduced models is a generalization of the so called Lax-Foldy reduced model.